program upp;

uses common,ppot;
Var
   i,j,k,l,o,m,n,dbg	 : integer;
   i1,i2,p1,p2,p3,ina	 : integer;
   txt			 : text;
   r			 : array[1..nsites] of v3;
   ene			 : array[1..nsites] of myreal;
   bl			 : array[1..3] of v3;
   itp			 : array[1..nsites] of integer;
   nat			 : integer;
   p			 : array[1..3] of integer;
   ar,u,ud,udd,xr,xe,ept : myreal;
   v,w			 : v3;
   b			 : boolean;
   ich			 : array[1..atps] of integer;
   ech			 : array[1..atps] of myreal;
   
function bound(z: myreal):myreal;
begin
   bound:=z;
   if (z>0.4999999) then bound:=z-1.0;
   if (z<-0.5000001) then bound:=z+1.0;
end; { bound }
begin
   dbg:=0;
   for i:=1 to 1000 do debug[i]:=false;
   pp_inp;
   for i:=1 to 3 do begin
      for o:=1 to 3 do read(bl[i,o]); readln;
   end;
   repeat
      nat:=nat+1;
      for o:=1 to 3 do read(r[nat,o]); readln(itp[nat]);
      if itp[nat]=0 then nat:=nat-1;
   until eof;
   epot:=0.0; ina:=0;
   for i:=1 to nat do ene[i]:=0.0;
   for i:=1 to nat do begin
      l:=0;
      for j:=i to nat do
	 for p1:=-iext[1] to iext[1] do
	    for p2:=-iext[2] to iext[2] do
	       for p3:=-iext[3] to iext[3] do begin
		  b:=((p1=0) and (p2=0) and (p3=0) and (i=j));
		  if not(b) then begin
		     for o:=1 to 3 do v[o]:=bound(r[j,o]-r[i,o]);
		     p[1]:=p1; p[2]:=p2; p[3]:=p3;
		     for m:=1 to 3 do w[m]:=0.0;
		     for m:=1 to 3 do
			for n:=1 to 3 do w[n]:=w[n]+(v[m]+p[m])*bl[m,n];
		     ar:=sqrt(w[1]*w[1]+w[2]*w[2]+w[3]*w[3]);
		     if (ar<rcut) then begin
			l:=l+1;
			upp(itp[i],itp[j],ar,u,ud,udd);
			epot:=epot+u;
			if (i=j) then epot:=epot-u/2.0;
			ene[i]:=ene[i]+u/2.0;
			if (i<>j) then ene[j]:=ene[j]+u/2.0;
		     end;
		  end;
	       end;
      if l>ina then ina:=l;
   end;
   for i:=1 to 110 do begin ich[i]:=0; ech[i]:=0.0; end;
   ept:=0.0;
   for i:=1 to nat do begin
      ich[itp[i]]:=ich[itp[i]]+1;
      ech[itp[i]]:=ech[itp[i]]+ene[i];
   end;
   writeln(' CHEM       N         x    E[ev/atom]');
   for i:=1 to 110 do if ich[i]>0 then begin
      xr:=ich[i]/nat;
      xe:=ech[i]/ich[i];
      writeln(i:5,ich[i]:8,xr:10:4,xe:10:4);
      ept:=ept+xe*xr;
   end;
   writeln(nat:1,' atoms,  E=',epot/nat:1:8,'[eV]');
   {writeln('Test E=',ept:1:8);}
   writeln('RCUT=',rcut:1:2,'A: max. ',ina:1,' atoms in the interaction sphere.');
end.
